Aufgabe ist es aus der Chemotion Repository (https://www.chemotion-repository.net/welcome) alle bisher hinterlegten Substanzen zu klassifizieren und eine deskriptive Statistik des erhaltenen Datensatzes anzufertigen. Hierzu wird im Folgendem das verwendete R-Script beschrieben.

Im ersten Schritt werden die benötigten Pakete geladen:

library(readr)
library(rinchi)
library(magrittr)
library(classyfireR)
library(dplyr)
library(sunburstR)
library(Rcpp)
library(RMassBank)
library(ggplot2)
library(plotly)
library(metfRag)

Die benötigten canonical smiles der Substanzen aus der Chemotion Repository können als XLSX-Datei heruntergeladen werden. Hierfür muss ein Account auf der Homepage der Chemotion Repository angelegt werden. Unter dem Reiter “My DB” können unter “Chemotion” alle bisher eingetragenen Substanzen markiert und dann als XLSX-Datei exportiert werden. In der heruntergeladenen XLSX-Datei befinden sich anschließend noch Duplikate, sowie Lösungsmittel (solvents) und Reaktionspartner (reactants). Nach Entfernung dieser Substanzen wird die XLSX-Datei in eine CSV-Datei konvertiert und kann dann in R eingelesen werden.

smiles <- read.csv("~/R/Chemotion/ChemotionViz/Data/sample_export_16.03.2021_8.12_noDup.csv")[ ,5]
head(smiles,3)
## [1] "NNc1ccc(cc1)Br.Cl"                                                             
## [2] "CC1=Nc2c(C1(C)C)cc(cc2)/C=C/C(C(C(C(C(C(C(F)(F)F)(F)F)(F)F)(F)F)(F)F)(F)F)(F)F"
## [3] "C[N+]1=C(C)C(c2c1ccc(c2)/C=C/C(C(C(C(F)(F)F)(F)F)(F)F)(F)F)(C)C.[I-]"

Mit dem Paket “rinchi” (https://rdrr.io/github/CDK-R/rinchi/man/getinchi.html#heading-0) werden nun die canonical smiles in InChikeys konvertiert, damit diese zur Klassifizierung notwenig sind.

inchikey <-sapply(smiles, get.inchi.key)

Im nächsten Schritt werden die InChiKeys verwendet, um eine Klassifizierung der Substanzen durchzuführen. Die Klassifizierung erfolgt mit dem Paket “classyfireR” (https://github.com/aberHRML/classyfireR). Die Funktion “get_classification” erzeugt eine Klassifizierungsliste. In dieser wird für jede Substanz ein Objekt welches Listen mit Tibbles enthält angelegt. Aus diesem Objekt können dann zum Beispiel Metadaten und die Klassifizierungstabelle ausgelesen werden (siehe folgende Ausgabe, Objekt 1). Für Substanzen die nicht klassifiziert werden konnten wird ein Objekt vom Typ “NULL” erzeugt (siehe folgende Ausgabe, Objekt 2+3).

Classification_List <- sapply(inchikey, get_classification)
head(Classification_List,3)
## $`NNc1ccc(cc1)Br.Cl`
## ── ClassyFire Object ───────────────────────────────────── classyfireR v0.3.6 ── 
## Object Size: 12.2 Kb 
##  
## Info: 
## ● InChIKey=RGGOWBBBHWTTRE-UHFFFAOYSA-N
##   
## ● Cl.NNC1=CC=C(Br)C=C1
##   
## ● Classification Version: 2.1
##   
## kingdom : Organic compounds
## └─superclass : Benzenoids
##   └─class : Benzene and substituted derivatives
##     └─subclass : Phenylhydrazines
## 
## $`CC1=Nc2c(C1(C)C)cc(cc2)/C=C/C(C(C(C(C(C(C(F)(F)F)(F)F)(F)F)(F)F)(F)F)(F)F)(F)F`
## NULL
## 
## $`C[N+]1=C(C)C(c2c1ccc(c2)/C=C/C(C(C(C(F)(F)F)(F)F)(F)F)(F)F)(C)C.[I-]`
## NULL

Nun kann der gewünschte Klassifizierungsgrad der Substanzen ausgewählt werden, hier level=“class” (“kingdom”=1; “superclass”=2; “class”=3; “subclass”=4).

level=3

Zum Auslesen des Klassifizierungsgrades wird eine For-Schleife genutzt. Für diese wird zunächst ein leerer Vektor mit der Länge der erzeugten Klassifizierungsliste erzeugt. Die For-Schleife iteriert über die einzelnen Objekte der Klassifizierungsliste, um den gewünschten Klassifizierungsgrad auszulesen. Mit einer If-Bedingung werden Objekte abgefangen, die vom Typ “NULL” sind. Die Schleife füllt den zuvor erzeugten leeren Vektor mit den Klassifizierungsgrad der jeweiligen Substanz.

class <- vector("character", length(Classification_List[]))
for (i in 1:length(Classification_List[])) {
  if( is.null(Classification_List[[i]]) == 'TRUE') { 
  } 
  else if (is.null(Classification_List[[i]]@classification[["Classification"]][level]) == 'TRUE'){
  }
  else { 
    class[i] <- Classification_List[[i]]@classification[["Classification"]][level]
  }
}

Anzahl der zu klassifizierenden Substanzen aus der Chemotion Repository (Größe des Datensatzes):

## [1] 2343

Im nächsten Schritt werden die leeren Strings aus dem class-Vektor in “NA” überführt und dann entfernt.

class[class==""] <- NA
class[!is.na(class)]

Anschließend wird der Klassifizierungsgrad der Substanzen alphabetisch sortiert, in einen Data Frame geschrieben und ausgezählt. Die Ausgabe ist dann ein Data Frame mit den alphabetisch sortierten Substanzklassen und der Anzahl dieser im Datensatz. Substanzklassen, die in ihrem Namen einen Bindestrich haben, werden im folgenden Sunburst-Plot nicht richtig geplotted, deshalb werden die Bindestriche in Unterstriche umgewandelt.

class_sorted <- sort(class)
class_dash <- gsub("-", "_", class_sorted)
df <- data.frame(class_dash)
(n_classes <- dplyr::count(df,class_dash))

Anzahl der erfolgreich klassifizierten Substanzen aus der Chemotion Repository:

## [1] 1060

Mithilfe dieses Data Frames kann die Darstellung des Klassifizierungslevel “class” aller klassifizierten Substanzen prozentual als Sunburst-Plot erfolgen. Die Erstellung des Sunburst Plots erfolgte mit dem Paket “sunburstR” (https://github.com/timelyportfolio/sunburstR)

sunburst(data = data.frame(n_classes), legend = FALSE)
Legend

Neben der Darstellung der prozentualen Anteile der Substanzklassen können mithilfe des Paketes “RMassBank” (https://bioconductor.org/packages/release/bioc/html/RMassBank.html) die molaren Massen der Substanzen aus den canonical smiles berechnet werden. Der erstellte Data Frame zeigt die ersten drei Substanzen aus dem Datensatz mit der jeweiligen exakten molaren Masse [g/mol] an.

substance_mass <- sapply(smiles, smiles2mass)
df_mass <- data.frame(substance_mass)
mass <- data.frame(substance= row.names(df_mass), df_mass, row.names=NULL)
head(mass,3)

Mit dem Paket “ggplot2” (https://www.rdocumentation.org/packages/ggplot2/versions/3.3.3) können die Substanzen mit der dazugehörigen molaren Masse dargestellt werden und mit dem Paket “ggplotly” (https://www.rdocumentation.org/packages/plotly/versions/4.9.3/topics/ggplotly) werden die Grafiken interaktiv dargestellt.

ggplot(mass,aes(substance,substance_mass))+
  geom_col(color="darkblue")+
  labs(title="Substances vs Molar mass",x="Substance", y = "Exact mass (g/mol)")

histogram <- ggplot(mass,aes(substance_mass))+
  geom_histogram(binwidth=20,color="darkblue", fill="lightblue")+
  labs(title="Histogram plot: Molar mass vs count",x="Exact mass (g/mol)", y = "Count")

ggplotly(histogram)

Beispiel-Plot eines hierarchischen Clusterings der Substanzen mit dem Paket ‘metfRag’:

mols <- parse.smiles(smiles)
dummy <- mapply(set.property, mols, "Score", c(1:2343))
plotCluster(mols, h=0.2)